This notebook looks at classifying tumor cells with
CellAssign in a Ewing sarcoma sample, SCPCS000491.
We use CellAssign with a variety of references:
We then compare the results from running CellAssign with
each reference to the tumor/normal classifications obtained from looking
only at marker genes. We also look at marker gene expression in the
assigned cells.
suppressPackageStartupMessages({
# load required packages
library(SingleCellExperiment)
library(ggplot2)
})
## Warning: replacing previous import 'S4Arrays::makeNindexFromArrayViewport' by
## 'DelayedArray::makeNindexFromArrayViewport' when loading 'SummarizedExperiment'
# The base path for the OpenScPCA repository, found by its (hidden) .git directory
repository_base <- rprojroot::find_root(rprojroot::is_git_root)
# The current data directory, found within the repository base directory
data_dir <- file.path(repository_base, "data", "current")
sample_dir <- file.path(data_dir, "SCPCP000015", params$sample_id)
# The path to this module
module_base <- file.path(repository_base, "analyses", "cell-type-ewings")
# source in helper functions make_jaccard_matrix() and jaccard()
jaccard_functions <- file.path(module_base, "scripts", "utils", "jaccard-functions.R")
source(jaccard_functions)
# Input files
sce_filename <- glue::glue("{params$library_id}_processed.rds")
sce_file <- file.path(sample_dir, sce_filename)
tumor_markers_file <- file.path(module_base, "references", "tumor-marker-genes.tsv")
visser_markers_file <- file.path(module_base, "references", "visser-all-marker-genes.tsv")
# output annotations
cellassign_annotations_file <- file.path(params$results_dir, glue::glue("{params$library_id}_cellassign-classifications.tsv"))
Define some helper functions for creating a reference matrix and
obtaining cell type annotations from CellAssign
predictions.
# function to get assigned cell type based on max prediction returned by CellAssign
get_celltype_assignment <- function(predictions){
# if no predictions, assign every cell to NA
if(length(colnames(predictions)) == 1){
celltype_assignments <- data.frame(
barcode = predictions$barcode,
celltype = NA
)
} else {
# get individual cell type assignments
# those with the max prediction
celltype_assignments <- predictions |>
tidyr::pivot_longer(
!barcode,
names_to = "celltype",
values_to = "score"
) |>
dplyr::group_by(barcode) |>
dplyr::slice_max(score, n = 1) |>
dplyr::ungroup()
}
return(celltype_assignments)
}
# read in processed sce
sce <- readr::read_rds(sce_file)
# read in tumor marker genes table
tumor_markers_df <- readr::read_tsv(tumor_markers_file) |>
# account for genes being from multiple sources
dplyr::select(cell_type, ensembl_gene_id, gene_symbol) |>
dplyr::distinct()
## Rows: 13 Columns: 4
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): cell_type, gene_symbol, ensembl_gene_id, source
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# read in all marker genes
visser_markers_df <- readr::read_tsv(visser_markers_file) |>
# account for genes being from multiple sources
dplyr::select(cell_type, ensembl_gene_id, gene_symbol) |>
dplyr::distinct()
## Rows: 18 Columns: 4
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): cell_type, gene_symbol, ensembl_gene_id, source
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# read in tumor normal classifications
# obtained from manual classification in `01-marker-gene-tumor-classifications.Rmd`
classifications_df <- readr::read_tsv(params$marker_gene_classification)
## Rows: 12749 Columns: 2
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (2): barcodes, marker_gene_classification
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
tumor_predictions <- readr::read_tsv(params$tumor_marker_predictions)
## Rows: 12749 Columns: 3
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): barcode
## dbl (2): tumor, other
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
visser_predictions <- readr::read_tsv(params$visser_marker_predictions)
## Rows: 12749 Columns: 6
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): barcode
## dbl (5): tumor, immune, mesenchymal-like cells, endothelial cells, other
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
panglao_predictions <- readr::read_tsv(params$panglao_predictions)
## Rows: 12749 Columns: 4
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): barcode
## dbl (3): Endothelial cells, Fibroblasts, other
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# check if marker gene annotations are present
if(all(is.na(classifications_df$marker_gene_classification))){
has_marker_gene <- FALSE
message("No annotations were available using only marker gene expression.
Any plots comparing CellAssign to marker gene annotation will be skipped.")
} else {
has_marker_gene <- TRUE
}
First let’s see if we can use CellAssign to classify
cells as tumor or normal using just the list of tumor marker genes.
We will read in the predictions file and find the predicted cell type for each cell, either “tumor” or “other”. Then we’ll show a heatmap of the scores for each cell and cell type. Included in that heatmap will be annotation of which cells were classified as tumor or normal using the manual classification.
# get cell type assignments
tumor_celltype_assignments <- get_celltype_assignment(tumor_predictions)
if(all(is.na(tumor_celltype_assignments$celltype))){
has_tumor_assignments <- FALSE
message("No CellAssign assignments found for tumor marker genes.")
} else {
# create a table of assignments
table(tumor_celltype_assignments$celltype)
has_tumor_assignments <- TRUE
}
# define annotations to use in all heatmaps
if(has_marker_gene){
# get annotations for heatmaps
annotation_df <- data.frame(
marker_gene_classification = classifications_df$marker_gene_classification,
row.names = classifications_df$barcodes
)
} else {
annotation_df <- NULL
}
# heatmap of prediction scores
tumor_predictions |>
tibble::column_to_rownames("barcode") |>
as.matrix() |>
pheatmap::pheatmap(show_rownames = FALSE,
annotation_row = annotation_df)
Let’s also look at the cell annotations on a UMAP and then look at the confusion matrix.
# rename cell type column
tumor_celltype_assignments <- tumor_celltype_assignments |>
dplyr::rename(cellassign_annotations = celltype)
# create a dataframe to use for plotting with UMAP and cell information
celltype_df <- sce |>
scuttle::makePerCellDF(use.dimred = "UMAP") |>
# replace UMAP.1 with UMAP1
dplyr::rename_with(
\(x) stringr::str_replace(x, "^UMAP\\.", "UMAP")
) |>
# join with tumor normal
dplyr::left_join(classifications_df, by = "barcodes") |>
# join with filtered cell type assignments from cellassign
dplyr::left_join(tumor_celltype_assignments, by = c("barcodes" = "barcode"))
ggplot(celltype_df, aes(x = UMAP1, y = UMAP2, color = cellassign_annotations)) +
geom_point(size = 0.5, alpha = 0.5) +
theme_bw()
# prep annotations for confusion matrix
celltype_df <- celltype_df |>
dplyr::mutate(cellassign_annotations = ifelse(
# make sure annotations match marker gene annotations
cellassign_annotations == "tumor", "Tumor", "Normal"
)) |>
# make tumor the positive class
dplyr::mutate(
cellassign_annotations = forcats::fct_relevel(cellassign_annotations, "Tumor"),
marker_gene_classification = forcats::fct_relevel(marker_gene_classification, "Tumor")
)
caret::confusionMatrix(
table(
celltype_df$marker_gene_classification,
celltype_df$cellassign_annotations)
)
## Confusion Matrix and Statistics
##
##
## Tumor Normal
## Tumor 1 6012
## Normal 0 6736
##
## Accuracy : 0.5284
## 95% CI : (0.5197, 0.5371)
## No Information Rate : 0.9999
## P-Value [Acc > NIR] : 1
##
## Kappa : 2e-04
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 1.000e+00
## Specificity : 5.284e-01
## Pos Pred Value : 1.663e-04
## Neg Pred Value : 1.000e+00
## Prevalence : 7.844e-05
## Detection Rate : 7.844e-05
## Detection Prevalence : 4.716e-01
## Balanced Accuracy : 7.642e-01
##
## 'Positive' Class : Tumor
##
To confirm that the cells with tumor cells have higher expression of
marker genes, let’s look at the mean expression of marker genes between
tumor and other cells as classified by CellAssign.
# get marker gene expression
marker_gene_exp <- logcounts(sce[tumor_markers_df$ensembl_gene_id, ]) |>
as.matrix() |>
t() |>
as.data.frame() |>
tibble::rownames_to_column("barcodes")
marker_gene_exp_df <- celltype_df |>
dplyr::left_join(marker_gene_exp, by = "barcodes") |>
tidyr::pivot_longer(
cols = starts_with("ENSG"),
names_to = "ensembl_gene_id",
values_to = "gene_expression"
) |>
dplyr::group_by(barcodes) |>
dplyr::mutate(mean_exp = mean(gene_expression))
ggplot(marker_gene_exp_df, aes(x = mean_exp, fill = cellassign_annotations)) +
geom_density() +
facet_grid(rows = vars(cellassign_annotations)) +
theme_bw()
We can also test to see if classifying tumor cells works with adding
a few more cell types, rather than just have the two options. Below we
will use CellAssign but with all tumor marker genes and all
marker genes in Visser et al. for non-tumor cells.
# get cell type assignments
visser_celltype_assignments <- get_celltype_assignment(visser_predictions)
if(all(is.na(visser_celltype_assignments$celltype))){
has_visser_assignments <- FALSE
message("No CellAssign assignments found for tumor marker genes.")
} else {
# create a table of assignments
table(visser_celltype_assignments$celltype)
has_visser_assignments <- TRUE
}
# heatmap of prediction scores
visser_predictions |>
tibble::column_to_rownames("barcode") |>
as.matrix() |>
pheatmap::pheatmap(show_rownames = FALSE,
annotation_row = annotation_df)
Below we will look at which cells are assigned to which cell type and
then compare assignments between CellAssign and manual
marker gene assignment.
# first rename column to combine with cell type df
visser_celltype_assignments <- visser_celltype_assignments |>
dplyr::rename(visser_cellassign_annotations = celltype)
celltype_df <- celltype_df |>
dplyr::left_join(visser_celltype_assignments, by = c("barcodes" = "barcode"))
# umap showing cell type annotations from using Visser reference
ggplot(celltype_df, aes(x = UMAP1, UMAP2, color = visser_cellassign_annotations)) +
geom_point(size = 0.5, alpha = 0.5) +
theme_bw()
# get jaccard similarity index
jaccard_mtx <- make_jaccard_matrix(celltype_df, "visser_cellassign_annotations", "marker_gene_classification")
ComplexHeatmap::Heatmap(
t(jaccard_mtx), # transpose because matrix rows are in common & we want a vertical arrangement
col = circlize::colorRamp2(c(0, 1), colors = c("white", "darkslateblue")),
border = TRUE,
## Row parameters
cluster_rows = FALSE,
row_title = "Marker gene",
row_title_gp = grid::gpar(fontsize = 12),
row_title_side = "left",
row_names_side = "left",
row_dend_side = "right",
row_names_gp = grid::gpar(fontsize = 10),
## Column parameters
cluster_columns = FALSE,
column_title = "CellAssign - Visser",
column_title_gp = grid::gpar(fontsize = 12),
column_names_side = "bottom",
column_names_gp = grid::gpar(fontsize = 10),
column_names_rot = 90,
## Legend parameters
heatmap_legend_param = list(
title = "Jaccard index",
direction = "vertical",
legend_width = unit(1.5, "in")
),
show_heatmap_legend = TRUE,
)
The one cell type that is identified here and in the
CellAssign annotations with the reference from
scpca-nf is Endothelial cells. Let’s compare the
classifications of Endothelial cells between the two references.
# first make sure cells have been classified as endothelial
has_endothelial <- "Endothelial cells" %in% unique(celltype_df$cellassign_celltype_annotation) &&
"endothelial cells" %in% unique(celltype_df$visser_cellassign_annotations)
# if present, compare across both references
if(has_endothelial){
# prep annotations for confusion matrix
celltype_df <- celltype_df |>
dplyr::mutate(
caret_visser_annotations = ifelse(
# make sure annotations match marker gene annotations
visser_cellassign_annotations == "endothelial cells", "Endothelial cells", "Other"
),
caret_cellassign_annotations = ifelse(
cellassign_celltype_annotation == "Endothelial cells", "Endothelial cells", "Other"
),
# make tumor the positive class
caret_visser_annotations = forcats::fct_relevel(caret_visser_annotations, "Endothelial cells"),
caret_cellassign_annotations = forcats::fct_relevel(caret_cellassign_annotations, "Endothelial cells")
)
caret::confusionMatrix(
table(
celltype_df$caret_cellassign_annotations,
celltype_df$caret_visser_annotations)
)
} else {
message("No endothelial cells were classified by `CellAssign`.")
}
## No endothelial cells were classified by `CellAssign`.
Below, we look at the individual expression of each marker genes in the cells. These plots are grouped to show all marker genes for a given cell type in a single panel.
# get all markers
all_markers <- visser_markers_df |>
dplyr::pull(ensembl_gene_id)
# get expression for all markers
marker_gene_exp <- logcounts(sce[all_markers, ]) |>
as.matrix() |>
t() |>
as.data.frame() |>
tibble::rownames_to_column("barcodes")
# add in gene expression to cell type df for plotting
marker_gene_exp_df <- celltype_df |>
dplyr::left_join(marker_gene_exp, by = "barcodes") |>
tidyr::pivot_longer(
cols = starts_with("ENSG"),
names_to = "ensembl_gene_id",
values_to = "gene_expression"
) |>
dplyr::left_join(visser_markers_df)
## Joining with `by = join_by(ensembl_gene_id)`
# get a list of all celltypes in the reference
celltypes <- unique(visser_markers_df$cell_type)
# for each cell type, plot each individual marker as a UMAP
celltypes |>
purrr::map(\(celltype){
plot_df <- marker_gene_exp_df |>
dplyr::filter(cell_type == celltype)
# faceted umap showing a umap panel for each marker gene
ggplot(plot_df, aes(x = UMAP1, y = UMAP2, color = gene_expression)) +
geom_point(alpha = 0.1, size = 0.2) +
facet_wrap(vars(gene_symbol)) +
scale_color_viridis_c() +
labs(
title = celltype,
color = "Log-normalized gene expression"
) +
# remove axis numbers and background grid
scale_x_continuous(labels = NULL, breaks = NULL) +
scale_y_continuous(labels = NULL, breaks = NULL) +
theme(
aspect.ratio = 1,
legend.position = "bottom",
axis.title = element_text(size = 9, color = "black"),
strip.text = element_text(size = 8),
legend.title = element_text(size = 9),
legend.text = element_text(size = 8)
) +
guides(colour = guide_colorbar(title.position = "bottom", title.hjust = 0.5)) +
theme_bw()
})
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
Below we look at the expression of marker genes for each cell type
across all cells identified to be each type by CellAssign.
We would expect to see expression of the marker genes for tumor cells to
be higher in the cells that are identified to be tumor cells compared to
other cell types.
grouped_celltype_df <- marker_gene_exp_df |>
# group by barcode and cell type to get the sum
dplyr::group_by(barcodes, cell_type) |>
dplyr::mutate(mean_exp = mean(gene_expression))
ggplot(grouped_celltype_df, aes(x = mean_exp, color = cell_type)) +
geom_density() +
facet_wrap(vars(visser_cellassign_annotations)) +
theme_bw() +
labs(
color = "Marker gene group"
)
The last thing we will do here is run CellAssign using some of the same reference that we used when running CellAssign as part of scpca-nf. In particular, the majority of non-muscle cells (which we believe to be tumor cells) were identified as either Endothelial cells or Fibroblasts. So we will use a reference that takes the markers for Endothelial cells and Fibroblasts from PanglaoDB and combines with the tumor markers.
# get cell type assignments
panglao_celltype_assignments <- get_celltype_assignment(panglao_predictions)
if(all(is.na(panglao_celltype_assignments$celltype))){
has_panglao_assignments <- FALSE
message("No CellAssign assignments found for tumor marker genes.")
} else {
# create a table of assignments
table(panglao_celltype_assignments$celltype)
has_panglao_assignments <- TRUE
}
# heatmap of prediction scores
panglao_predictions |>
tibble::column_to_rownames("barcode") |>
as.matrix() |>
pheatmap::pheatmap(show_rownames = FALSE,
annotation_row = annotation_df)
# first rename column to combine with cell type df
panglao_celltype_assignments <- panglao_celltype_assignments |>
dplyr::rename(panglao_cellassign_annotations = celltype)
celltype_df <- celltype_df |>
dplyr::left_join(panglao_celltype_assignments, by = c("barcodes" = "barcode"))
# umap showing cell type annotations from using Visser reference
ggplot(celltype_df, aes(x = UMAP1, UMAP2, color = panglao_cellassign_annotations)) +
geom_point(size = 0.5, alpha = 0.5) +
theme_bw()
Let’s directly compare these annotations to the
CellAssign annotations obtained from
scpca-nf.
# get jaccard similarity index
jaccard_mtx <- make_jaccard_matrix(celltype_df, "cellassign_celltype_annotation", "panglao_cellassign_annotations")
ComplexHeatmap::Heatmap(
t(jaccard_mtx), # transpose because matrix rows are in common & we want a vertical arrangement
col = circlize::colorRamp2(c(0, 1), colors = c("white", "darkslateblue")),
border = TRUE,
## Row parameters
cluster_rows = FALSE,
row_title = "Panglao + marker genes",
row_title_gp = grid::gpar(fontsize = 12),
row_title_side = "left",
row_names_side = "left",
row_dend_side = "right",
row_names_gp = grid::gpar(fontsize = 10),
## Column parameters
cluster_columns = FALSE,
column_title = "scpca-nf CellAssign",
column_title_gp = grid::gpar(fontsize = 12),
column_names_side = "bottom",
column_names_gp = grid::gpar(fontsize = 10),
column_names_rot = 90,
## Legend parameters
heatmap_legend_param = list(
title = "Jaccard index",
direction = "vertical",
legend_width = unit(1.5, "in")
),
show_heatmap_legend = TRUE,
)
celltype_df <- celltype_df |>
dplyr::select(barcodes,
marker_gene_cellassign_annotations = cellassign_annotations,
visser_cellassign_annotations,
panglao_cellassign_annotations)
readr::write_tsv(celltype_df, cellassign_annotations_file)
# record the versions of the packages used in this analysis and other environment information
sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices datasets utils methods
## [8] base
##
## other attached packages:
## [1] ggplot2_3.5.1 SingleCellExperiment_1.26.0
## [3] SummarizedExperiment_1.34.0 Biobase_2.64.0
## [5] GenomicRanges_1.56.0 GenomeInfoDb_1.40.0
## [7] IRanges_2.38.0 S4Vectors_0.42.0
## [9] BiocGenerics_0.50.0 MatrixGenerics_1.16.0
## [11] matrixStats_1.3.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 jsonlite_1.8.8
## [3] shape_1.4.6.1 magrittr_2.0.3
## [5] farver_2.1.2 rmarkdown_2.26
## [7] GlobalOptions_0.1.2 zlibbioc_1.50.0
## [9] vctrs_0.6.5 DelayedMatrixStats_1.26.0
## [11] htmltools_0.5.8.1 S4Arrays_1.4.0
## [13] forcats_1.0.0 SparseArray_1.4.3
## [15] pROC_1.18.5 caret_6.0-94
## [17] sass_0.4.9 parallelly_1.37.1
## [19] bslib_0.7.0 plyr_1.8.9
## [21] lubridate_1.9.3 cachem_1.0.8
## [23] lifecycle_1.0.4 iterators_1.0.14
## [25] pkgconfig_2.0.3 Matrix_1.7-0
## [27] R6_2.5.1 fastmap_1.2.0
## [29] GenomeInfoDbData_1.2.12 future_1.33.2
## [31] clue_0.3-65 digest_0.6.35
## [33] colorspace_2.1-0 rprojroot_2.0.4
## [35] beachmat_2.20.0 labeling_0.4.3
## [37] fansi_1.0.6 timechange_0.3.0
## [39] httr_1.4.7 abind_1.4-5
## [41] compiler_4.4.0 proxy_0.4-27
## [43] bit64_4.0.5 withr_3.0.0
## [45] doParallel_1.0.17 BiocParallel_1.38.0
## [47] highr_0.10 MASS_7.3-60.2
## [49] lava_1.8.0 DelayedArray_0.30.1
## [51] rjson_0.2.21 ModelMetrics_1.2.2.2
## [53] tools_4.4.0 future.apply_1.11.2
## [55] nnet_7.3-19 glue_1.7.0
## [57] nlme_3.1-164 grid_4.4.0
## [59] cluster_2.1.6 reshape2_1.4.4
## [61] generics_0.1.3 recipes_1.0.10
## [63] gtable_0.3.5 tzdb_0.4.0
## [65] class_7.3-22 tidyr_1.3.1
## [67] data.table_1.15.4 hms_1.1.3
## [69] utf8_1.2.4 XVector_0.44.0
## [71] foreach_1.5.2 pillar_1.9.0
## [73] stringr_1.5.1 vroom_1.6.5
## [75] circlize_0.4.16 splines_4.4.0
## [77] dplyr_1.1.4 lattice_0.22-6
## [79] renv_1.0.7 survival_3.5-8
## [81] bit_4.0.5 tidyselect_1.2.1
## [83] ComplexHeatmap_2.20.0 scuttle_1.14.0
## [85] knitr_1.46 xfun_0.44
## [87] hardhat_1.3.1 timeDate_4032.109
## [89] pheatmap_1.0.12 stringi_1.8.4
## [91] UCSC.utils_1.0.0 yaml_2.3.8
## [93] evaluate_0.23 codetools_0.2-20
## [95] tibble_3.2.1 BiocManager_1.30.23
## [97] cli_3.6.2 rpart_4.1.23
## [99] munsell_0.5.1 jquerylib_0.1.4
## [101] Rcpp_1.0.12 globals_0.16.3
## [103] png_0.1-8 parallel_4.4.0
## [105] gower_1.0.1 readr_2.1.5
## [107] sparseMatrixStats_1.16.0 listenv_0.9.1
## [109] viridisLite_0.4.2 ipred_0.9-14
## [111] scales_1.3.0 prodlim_2023.08.28
## [113] e1071_1.7-14 purrr_1.0.2
## [115] crayon_1.5.2 GetoptLong_1.0.5
## [117] rlang_1.1.3